#ifndef TEMPLATEFIT_H
#define TEMPLATEFIT_H

#ifndef MATHFIT_H
#include "MathFit.h"
#endif

class TemplateFit
{
 public:
  static void Fit(const TH1 *hdata0, const Int_t ntmpl,  TH1D * htmpls[], const TString tag, const Bool_t kfit, Double_t *pars);

 private:
  static Double_t TemplateFitFCN(const Double_t *xx, const Double_t *par);

  static void DrawFit(TH1D * hdata, const Double_t *par, const TString tag);

  static Int_t fNTmpl;
  static TH1D *fTmpls[10];
  static BFUNC fFunc;
};

Int_t TemplateFit::fNTmpl = 0;
TH1D * TemplateFit::fTmpls[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
BFUNC TemplateFit::fFunc = 0x0;

void TemplateFit::Fit(const TH1 *hdata0, const Int_t ntmpl,  TH1D * htmpls[], const TString tag, const Bool_t kfit, Double_t *pars)
{
  TH1D *hdata = (TH1D*) hdata0->Clone("templatefithdata0");

  if(ntmpl>10){
    printf("TemplateFit::Fit ntmpl>10! %d\n", ntmpl);exit(1);
  }

  fNTmpl =  ntmpl;
  for(Int_t ii=0; ii<fNTmpl; ii++){
    fTmpls[ii] = (TH1D*) htmpls[ii]->Clone(Form("templatefithtemp%d",ii));
  }

  fFunc = TemplateFitFCN;
  if(kfit){
    MathFit::InputData(hdata);

    const Int_t npar = fNTmpl;
    Double_t errs[10];

    const TString report = MathFit::MnFit(fFunc, npar, pars, errs);
  }

  DrawFit(hdata, pars, tag);

  delete hdata;
  for(Int_t ii=0; ii<fNTmpl;ii++){
    delete fTmpls[ii];
    fTmpls[ii]=0x0;
  }
}

void TemplateFit::DrawFit(TH1D * hdata, const Double_t *pars, const TString tag)
{
  //hdata drawn as original

  TCanvas *c1=new TCanvas;
  gStyle->SetOptTitle(1);

  //hfit only uses binning info of hdata
  //hfit content set by fFunc according to pars and fHistFit
  TH1D *hfit = (TH1D*)hdata->Clone("templatefitdraw");
  for(Int_t ii=0; ii<=hfit->GetNbinsX()+1; ii++){
    const Double_t xi = hfit->GetBinCenter(ii);
    const Double_t yi = fFunc(&xi, pars);
    hfit->SetBinContent(ii, yi);
  }

  hdata->SetMarkerStyle(22);
  hdata->SetMarkerSize(1);
  hfit->SetLineColor(kRed);

  const TString tit(Form("points: %s, line fit", hdata->GetTitle()));
  hdata->SetTitle(tit);

  gPad->SetLogy(0);
  hdata->Draw("e");
  hfit->Draw("hist same");
  c1->Print(tag+"logy0.png");

  gPad->SetLogy();
  hdata->Draw("e");
  hfit->Draw("hist same");
  c1->Print(tag+"logy1.png");

  delete hfit;
  delete c1;
}

Double_t TemplateFit::TemplateFitFCN(const Double_t *xx, const Double_t *par)
{
  Double_t sumy = 0.0;
  for(Int_t ii=0; ii<fNTmpl; ii++){
    const Int_t ib = fTmpls[ii]->FindBin(xx[0]);
    sumy += par[ii]*(fTmpls[ii]->GetBinContent(ib));
  }
  return sumy;
}

#endif
